/* ParticleSet.h
 *
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#if !defined(_ParticleSet_)
#define _ParticleSet_

#include <boost/operators.hpp>
#include <vector>
#include <src/miscellaneous/OS_basic.h>
#include <src/data/Matrix.h>
/*! \file ParticleSet.h
    \brief One and two dimensional array with full algebra set of operations.

	These classes are conserved for compatibility and test purposes for the classes OS_Array.
	OS_Array classes are faster.
*/




template <class T, class P>
struct particle_helper
{
	typedef typename T particles_list;
	typedef typename P particle_type;
	/*! \brief check structure for particle
		particle copy constructor uses memcpy to speed up copy by value.
		This class check whether the base type supports memcpy.
		A compile time error occurs if particle::base_value_type (T in particle<T>) is not of fundamental type.
	*/
	template<class P, class B>
	struct check_fundamental {
		OS_STATIC_CHECK(false);
	};
	template<class P>
	struct check_fundamental<P,OS_true>
	{};
	typedef typename boost::is_fundamental<P>::type bool_type;
	typedef check_fundamental<P,bool_type> check;
};
/*! \brief one dimensional vector with algebra operations.
*/
template<class T>
class particle : public particle_helper<particle<T>, T>,
	private      
		boost::additive< particle<T>          // +,-
	,	boost::multipliable < particle<T>, T //*, /
	,	boost::bitwise<particle<T>              // ^(inner scalar product), !(undefined at present time)
 	,	boost::shiftable<particle<T>, std::size_t     // <<, >>
     > > > >
{
public :
	enum :std::size_t {Dim = 1};
		typedef particle					type;
		typedef Protected					Boundary_type;
		typedef OS_int						indice_type;
		typedef OS_size						size_type;
		typedef std::vector<T>				Base;
		typedef T *							iterator;
		typedef const T *					const_iterator;
		typedef T							value_type;
		typedef T							base_value_type;
		typedef T &							reference;
		typedef const T&					const_reference;

		virtual ~particle (){
			if (_Myfirst) delete _Myfirst;
		}
		particle() : _Myfirst(0), _size(0) {};

		explicit particle(size_type n, T* value = 0) : _size(n) {
			if (_size == 0) 
				_Myfirst = 0;
			else {
				_Myfirst = new T[n];
				if (value) {
					for (indice_type i=0;i<(indice_type)_size;++i)
						(*this)[i] = value[i];
				}
			}
		}
		particle(const particle &from)
		{
				_size = from.size();
			OS_DYNAMIC_CHECK( _size >= 0, "negative size detected in particle allocator");

			if (_size == 0) {
				_Myfirst = 0;
			}
			else
			{
				_Myfirst = new T[_size];
				memcpy(_Myfirst,from._Myfirst,_size*sizeof(T));
			}
		};
		particle(const std::vector<T > &from)
		{
			_size = from.size();
			_Myfirst = new T[_size];
			for (size_type i=0;i<_size;++i, ++_Myfirst)
				(*_Myfirst) = from[i];
		};
		particle& operator = (const std::vector<T > &from) {
			if ( _size != from.size() )
				resize(_size);
			for (size_type i=0;i<_size;++i, ++_Myfirst)
				(*_Myfirst) = from[i];
			return *this;
         };
		inline size_type dim() const {return 1;};
		inline size_type size() const {
			return _size;
		};
		inline size_type size_list() const {
			return _size;
		};


		inline T* begin() {
			return _Myfirst;
		}
		inline T* end() {
			return _Myfirst + _size;
		}
		inline const_iterator begin() const {
			return _Myfirst;
		}
		inline const_iterator end() const {
			return _Myfirst + _size;
		}
		inline T & operator[](indice_type n) {

			OS_DYNAMIC_CHECK(n>=0 && n<(indice_type)_size, "indice out of range in particle accessor");

			return *(_Myfirst + n);
		}
		inline T & operator[](indice_type n) const {

			OS_DYNAMIC_CHECK(n>=0 && n<(indice_type)_size, "indice out of range in particle accessor");

			return *(_Myfirst + n);
		}
		inline void assign (T from) {

			T* Ptr_ = _Myfirst;
			T* End_ = _Myfirst+_size;
  			for (; Ptr_< End_; ++Ptr_)
					(*Ptr_) = from;
            };
		inline particle& operator = (T from) {
			T* Ptr_ = _Myfirst;
			T* End_ = _Myfirst+_size;
  			for (; Ptr_< End_; ++Ptr_)
					(*Ptr_) = from;
			return *this;
            };

		inline type& operator <<= (std::size_t how) {
			if (size() > 1)
				memcpy(_Myfirst, _Myfirst+1,(size()-1)*sizeof(base_value_type));
			(*this)[size()-1]= OS_MAX_REAL;
			return *this;
		};
		inline type& operator >>= (std::size_t how) {
			if (size() > 1)
				memcpy(_Myfirst+1, _Myfirst,(size()-1)*sizeof(base_value_type));
			(*this)[0]= OS_MAX_REAL;
			return *this;
		};


		inline void resize(size_type n)
		{

			OS_DYNAMIC_CHECK(n>=0,"negative size inputed in particle resize");

			if (_Myfirst) delete _Myfirst;
			_size = n;
			if (_size == 0) 
				_Myfirst = 0;
			else
				_Myfirst = new T[n];
		}
		inline particle operator - () {
			particle return_(*this);
			T* Ptr_ = return_._Myfirst;
			T* End_ = Ptr_ +_size;
  			for (; Ptr_< End_; ++Ptr_)
				(*Ptr_) *= -1;
			return return_;
            };
		inline particle operator ! () {
			particle return_(*this);
			T* Ptr_ = return_._Myfirst;
			T* End_ = Ptr_ +_size;
  			for (; Ptr_< End_; ++Ptr_)
				(*Ptr_) = 1./(*Ptr_);
			return return_;
            };

		inline particle& operator = (const particle &from) {
		    if (_size != from.size())
				resize(from.size());
			memcpy(_Myfirst,from._Myfirst,_size*sizeof(T));
			return *this;
            };

       inline particle& operator+=(const particle &from)  {

			OS_DYNAMIC_CHECK(_size == from.size(), "sizes not comparable in particle inner addition"); 

			T* Ptr_ = _Myfirst;
			T* Ptr_From = from._Myfirst;
			T* End_ = _Myfirst+_size;
			T* End_From = Ptr_From +from._size;
  			for (; Ptr_ < End_ && Ptr_From < End_From; ++Ptr_, ++Ptr_From)
					(*Ptr_) += (*Ptr_From);
			return *this;
        };
       inline particle& operator-=(const particle &from)  {

			OS_DYNAMIC_CHECK(_size == from.size(), "sizes not comparable in particle inner substraction");

			T* Ptr_ = _Myfirst;
			T* Ptr_From = from._Myfirst;
			T* End_ = _Myfirst+_size;
			T* End_From = Ptr_From +from._size;
  			for (; Ptr_< End_ && Ptr_From < End_From; ++Ptr_, ++Ptr_From)
					(*Ptr_) -= (*Ptr_From);
			return *this;
        };
       inline particle& operator^=(const particle &from)  {

			OS_DYNAMIC_CHECK(_size == from.size(), "sizes not comparable in particle inner addition"); 

			T* Ptr_ = _Myfirst;
			T* Ptr_From = from._Myfirst;
			T* End_ = _Myfirst+_size;
			T* End_From = Ptr_From +from._size;
  			for (; Ptr_ < End_ && Ptr_From < End_From; ++Ptr_, ++Ptr_From)
					(*Ptr_) *= (*Ptr_From);
			return *this;
        };
       inline bool operator ==(const particle &from)  {
			if (_size != from._size) return false;
			T* Ptr_ = _Myfirst;
			T* Ptr_From = from._Myfirst;
			T* End_ = _Myfirst+_size;
			for (; Ptr_< End_ ; ++Ptr_, ++Ptr_From) {
				if ( (*Ptr_) != (*Ptr_From)) return false;
			};
			return true;
        };
	   inline bool operator !=(const particle &from)  { 
			return  !(*this == from);
	   };

        inline particle& operator*=(const T from)  {
			T* Ptr_ = _Myfirst;
			T* End_ = _Myfirst+_size;
  			for (; Ptr_< End_; ++Ptr_)
					(*Ptr_) *= from;
			return *this;
         };
        inline particle& operator/=(const T from)  {
			T* Ptr_ = _Myfirst;
			T* End_ = _Myfirst+_size;
  			for (; Ptr_< End_; ++Ptr_)
					(*Ptr_) /= from;
			return *this;
         };
        particle sup_bound(const particle &from)  {
			particle return_(size());
			T* Ptr_ = _Myfirst;
			T* Ptr_From = from._Myfirst;
			T* Ptr_return = return_._Myfirst;
			T* End_ = _Myfirst+_size;
  			for (; Ptr_< End_; ++Ptr_, ++Ptr_From,++Ptr_return )
				(*Ptr_return) = (*Ptr_From) > (*Ptr_) ? (*Ptr_From) : (*Ptr_);
			return return_;
         };
        particle inf_bound(const particle &from)  {
			particle return_(_size);
			T* Ptr_ = _Myfirst;
			T* Ptr_From = from._Myfirst;
			T* Ptr_return = return_._Myfirst;
			T* End_ = _Myfirst+_size;
  			for (; Ptr_< End_; ++Ptr_, ++Ptr_From,++Ptr_return )
				(*Ptr_return) = (*Ptr_From) > (*Ptr_) ? (*Ptr_) : (*Ptr_From);
			return return_;
         };
        particle projector(size_type dimension)  const {
			particle return_(dimension);
			return_.assign(0.);
			T* Ptr_ = _Myfirst;
			T* Ptr_return = return_._Myfirst;
			T* End_ = _Myfirst+std::min(dimension, _size);
  			for (; Ptr_< End_; ++Ptr_, ++Ptr_return )
				(*Ptr_return) = (*Ptr_);
			return return_;
         };
		static std::string get_name() {return "-1D-Particle-";};

		T* _Myfirst;
		size_type _size;
};


/*! \brief list of one dimensional vector with algebra operations.
*/
template<class T >
class particles_list : public std::vector<T>,
	private      
		boost::additive< particles_list<T>          // +,-
	,	boost::multipliable < particles_list<T>, typename T::base_value_type  //*
	,	boost::bitwise<particles_list<T>              // ^(inner scalar product), !(undefined at present time)
      > > > //> 
{

public :
		enum :std::size_t {Dim = 2};
		typedef typename T::base_value_type							base_value_type;
		typedef Protected											Boundary_type;
		typedef std::vector<T>										Base;
		typedef typename Base::value_type							particle_type;
		typedef typename Base::iterator								comp_iterator;
		typedef typename T::iterator								part_iterator;
		typedef OS_size												size_type;
		typedef OS_int												indice_type;
		typedef std::pair<OS_int,OS_int>							indice_list_type;
		typedef T&													reference;
		typedef const T&											const_reference;

		particles_list (const T &from) {
			my_resize(from.size(),1);
	   			for (size_type i = 0; i<size();++i)
					(*this)[i][0] = from[i];

		}


		particles_list & assign(base_value_type double_value) {
   			for (size_type i = 0; i<size();++i)
			{
					(*this)[i].assign(double_value);
			}
			return *this;
		}
		particles_list & operator = (base_value_type double_value) {
   			for (size_type i = 0; i<size();++i)
			{
					(*this)[i] = double_value;
			}
			return *this;
		}
		void my_resize(size_type nb, size_type dim) {
			particles_number_=nb;
			particles_dim=dim;
			resize(particles_number_);
			comp_iterator begin1 = begin(), end1 = end();
			for (;begin1 != end1;++begin1) { 
				begin1->resize(particles_dim);};
		};


   		particles_list(const Base &from)
		{
			my_resize(from.size(),from.begin()->size());
   			for (size_type i = 0; i<size();++i)
			{
				OS_DYNAMIC_CHECK(from[i].size() == particles_dim, "trying to build a list of particles from a vector of non constant particles size");
	   			for (size_type j = 0; j<particles_dim;++j)
					(*this)[i][j] = from[i][j];
			}
		};


		particles_list(size_type number_of_particles, size_type dim_of_particles)	: particles_number_(number_of_particles), particles_dim(dim_of_particles), Base(number_of_particles) 
		{
			comp_iterator begin1 = begin(), end1 = end();
			for (;begin1 != end1;++begin1) { 
				begin1->resize(particles_dim);};
		}
		particles_list(const indice_list_type from)	: particles_number_(from.first), particles_dim(from.second), Base(from.first) 
		{
			comp_iterator begin1 = begin(), end1 = end();
			for (;begin1 != end1;++begin1) { 
				begin1->resize(particles_dim);};
		}
		particles_list() : particles_number_(0), particles_dim(0) {};

		void push_back(const particle_type& p) {
			++particles_number_;
			if (dim() == 0)
				particles_dim = p.size();
			else
				OS_DYNAMIC_CHECK(p.size() == dim(), "incoherent");
			Base::push_back(p);
		};
		// to optimize using iterators. There is also an index incoherency with OS_Vector to study
		template <class particle_like>
		particles_list& operator = (const particle_like &from) {
		    my_resize(from.size(1),from.size(2));
   			for (size_type j = 0; j<size(2);++j)
   				for (size_type i = 0; i<size(1);++i)
					(*this)[j][i] = from[i][j];
			return *this;
            };
		particles_list& operator = (const particles_list &from) {
		    my_resize(from.size(),from.dim());
   			for (size_type i = 0; i<size();++i)
					(*this)[i] = from[i];
			return *this;
            };
		inline particles_list& operator+=(const particles_list &from)  {

			OS_DYNAMIC_CHECK(size() == from.size(), "sizes not comparable in particles_list inner addition" );

			for (size_type i = 0; i<size();++i)
					(*this)[i] += from[i];
			return *this;
		};
		inline particles_list& operator^=(const particles_list &from)  {
					OS_DYNAMIC_CHECK(from1.size() == from2.size() && from1.dim() == from2.dim(), "exterior product of inconsistent list of particles");
					particles_list result(from1.size(),1);
   					for (size_type i = 0; i< from1.size() ;++i)
						result[i] = (*this)[i]*from2[i];
					*this = result;
					return *this;
		};

		inline bool operator ==(const particles_list &from)  {
		   if (size() != from.size()) return false;
		   for (size_type i = 0; i<size();++i) {
					if ((*this)[i] != from[i]) return false;
		   };
		   return true;
		};

       inline particles_list& operator *=(const T &from)  {
			OS_DYNAMIC_CHECK(size() == from.size(), "sizes not comparable in particles_list outer multiplication" );
   			for (size_type i = 0; i<size();++i)
					(*this)[i] *= from[i];
			return *this;
        };
		inline particles_list operator -()  const {
			particles_list result(*this);
			iterator begin = result.begin();
			iterator end = result.end();
			for (;begin != end; ++begin)
				*begin *= -1.;
			return result;
		 };

       inline particles_list& operator /=(const T &from)  {

			OS_DYNAMIC_CHECK(size() == from.size(), "sizes not comparable in particles_list outer division" );

   			for (size_type i = 0; i<size();++i)
					(*this)[i] /= from[i];
			return *this;
        };


       inline particles_list& operator-=(const particles_list &from)  {

			OS_DYNAMIC_CHECK(size() == from.size(), "sizes not comparable in particles_list inner substraction" );

   			for (size_type i = 0; i<size();++i)
					(*this)[i] -= from[i];
			return *this;
        };
        inline particles_list& operator*=(const base_value_type from)  {
   			for (size_type i = 0; i<size();++i)
					(*this)[i] *= from;
			return *this;
         };
        inline particles_list& operator/=(const base_value_type from)  {
   			for (size_type i = 0; i<size();++i)
					(*this)[i] /= from;
			return *this;
         };
		particles_list projector(size_type dimension) const {
			if (dimension == dim())
				return *this;
			particles_list return_(indice_list_type(size(),dimension));
   			for (size_type i = 0; i< size() ;++i)
				return_[i] = (*this)[i].projector(dimension);
			return return_;
		}
        T sup_bound() const {
			T return_(dim());
   			for (size_type i = 0; i< dim() ;++i)
				return_[i] = -1e+99;
   			for (size_type i = 0; i< size() ;++i)
				return_ = return_.sup_bound((*this)[i]);
			return return_;
         };
        T inf_bound() const {
			T return_(dim());
   			for (size_type i = 0; i< dim() ;++i)
				return_[i] = 1e+99;
   			for (size_type i = 0; i< size() ;++i)
				return_ = return_.inf_bound((*this)[i]);
			return return_;
         };


		size_type size() const {return particles_number_;} ;
		indice_list_type size_list() const {return indice_list_type(size(),dim());}; 
		size_type size(OS_size i) const {
			if (i==2)
				return size();
			if (i==1)
				return dim();
			OS_DYNAMIC_CHECK(false,"particles_list : this is a two dimensional array type");
			return 0;
		} ;
		size_type dim() const {return particles_dim;};
		size_type begin_indices() const {return 0;}
		size_type end_indices() const {return particles_number_;};
		static std::string get_name() {return "-2D-Particle List-";};

private :
		size_t particles_number_;
		size_t particles_dim;
};

/*! \brief scalar product of one dimensional vector
*/
template<class T>
inline T operator *(const particle<T> & from1, const particle<T> & from2) {
			T product = 0;

			OS_DYNAMIC_CHECK(from1.size() == from2.size(),"scalar product of inconsistent particles");
	
			particle<T>::const_iterator Ptr_Left = from1.begin();
			particle<T>::const_iterator Ptr_Right = from2.begin();
			particle<T>::const_iterator End_ = from1.end();
  			for (; Ptr_Left != End_; ++Ptr_Left, ++Ptr_Right )
				product += (*Ptr_Left)*(*Ptr_Right);
			return product;
};
/*! \brief division of of one dimensional vector by a scalar
*/
template<class T>
inline particle<T> operator /(const particle<T> & from1, const typename particle<T>::base_value_type & from2) {
		typedef particle<T> A;
		A temp(from1);
		temp /= from2;
		return temp;
};
/*! \brief division of of one dimensional vector by a scalar
*/
template<class T>
inline particle<T> operator /(const typename particle<T>::base_value_type & from2, const particle<T> & from1) {
		typedef particle<T> A;
		typedef A::iterator iterator;
		A temp(from1);
		iterator beg = temp.begin();
		iterator end = temp.end();
		for (;beg != end; ++beg)
			*beg = from2 / (*beg);
		return temp;
};

/*! \brief scalar product of one dimensional vectors
*/
template<class T>
inline typename T::base_value_type operator *(const particles_list<T> & from1, const particles_list<T> & from2) {
			T::base_value_type product = 0;
			typedef particles_list<T>::size_type size_type;
			OS_DYNAMIC_CHECK(from1.size() == from2.size() && from1.dim() == from2.dim(), "scalar product of inconsistent list of particles");
   			for (size_type i = 0; i< from1.size() ;++i)
				product += from1[i]*from2[i];
			return product;
};


/*! \brief division of two dimensional vector by a scalar
*/
template<class T>
inline particles_list<T> operator /(const particles_list<T> & from1, const typename particles_list<T>::base_value_type & from2) {
		typedef particles_list<T> A;
		A temp(from1);
		temp /= from2;
		return temp;
};
/*! \brief division of two dimensional vector by a scalar
*/
template<class T>
inline particles_list<T> operator /(const typename particles_list<T>::base_value_type & from2, const particles_list<T> & from1) {
		typedef particles_list<T> A;
		A temp(from1);
		temp *= from2;
		return 1./temp;
};



/*! \brief boolean binary operator < for one dimensional vectors
*/
template<class T>
inline bool operator < (const particle<T> & from1, const particle<T> & from2) {

			OS_DYNAMIC_CHECK(from1.size() == from2.size(),"particle: size not comparable during comparaison");
	
			particle<T>::const_iterator Ptr_Left = from1.begin();
			particle<T>::const_iterator Ptr_Right = from2.begin();
			particle<T>::const_iterator End_ = from1.end();
  			for (; Ptr_Left != End_; ++Ptr_Left, ++Ptr_Right )
				if ( (*Ptr_Left) >= (*Ptr_Right))
					return false;
			return true;
};

/*! \brief power of two dimensional vectors
*/
template<typename T, typename A>
inline particles_list<T> pow(const particles_list<T> & from1, const A& from2) {
	pas encore teste.
	typedef particles_list<T>::iterator iterator;
	typedef particles_list<T>::const_iterator const_iterator;
	particles_list<T> temp(from1.size(1),from1.size(2));
	const_iterator beg = from1.begin();const_iterator end = from1.end();iterator temp_beg = temp.begin();
	for (;beg != end;beg++,temp_beg++)
		*temp_beg = pow(*beg,from2);
	return temp;
};
/*! \brief power of one dimensional vectors
*/
template<typename T, typename A>
inline particle<T> pow(const particle<T> & from1, const A& from2) {
	typedef particle<T>::iterator iterator;
	typedef particle<T>::const_iterator const_iterator;
	particle<T> temp(from1.size());
	const_iterator beg = from1.begin();const_iterator end = from1.end();iterator temp_beg = temp.begin();
	for (;beg != end;beg++,temp_beg++)
		*temp_beg = pow(*beg,from2);
	return temp;
};



/*! \brief serialization of one dimensional vectors
*/
template<class T>
std::istream &operator>>(std::istream &in, particle<T> &to) {
	char c;
	in >> c;
	std::vector<T> t;
  if (c != '(') {
    in.setstate(std::ios_base::failbit);
    return in;
  }
   do {
		  T temp;
          in >> temp;
		  in >> c;
		  t.push_back(temp);
		} while (c == ',');
  if (c != ')') in.setstate(std::ios_base::failbit);
  to = t;
  return in;
};

/*! \brief serialization of two dimensional vectors
*/
template<class T>
std::istream &operator>>(std::istream &in, particles_list<T> & t)
{
	particles_list<T>::Base temp;
	in >> temp;
	t = temp;
	return in;
};

/*! \brief serialization of one dimensional vectors
*/
template<class T>
std::ostream &operator<<(std::ostream &o, const particle<T> &from) {
        o << "[ ";
        for (OS_size j = 0; j < from.size(); j ++) {
              o << from[j];
			 if (j!= from.size()-1)
				 o << ",";
        };
        o << " ]";
        return o;
};
/*! \brief serialization of one dimensional vectors
*/
template <class T>
std::ostream &operator<<(std::ostream &o, const particles_list<T> &from) {
        o << "[ ";
        for (OS_size j = 0; j < from.size(); j ++) {
              o << from[j];
			 if (j!= from.size()-1)
				 o << ",";
        };
        o << " ]";
        return o;
};


#endif 
